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Summary 

Several paradigms for accessing Computer Aided Design (CAD) 
geometry during surface meshing for Computational Fluid Dynam- 
ics (CFD) are discussed. File translation, inconsistent geometry 
engines and non-native point construction are all identified as 
sources of non-robustness. The paper argues in favor of accessing 
CAD parts and assemblies in their native format, without transla- 
tion, and for the use of CAD-native predicates and constructors in 
surface mesh generation. The discussion also emphasizes the impor- 
tance of examining the computational requirements for exact evalu- 
ation of triangulation predicates during surface meshing. 

The native approach is demonstrated through an algorithm for the 
generation of closed manifold surface triangulations from CAD 
geometry. CAD parts and assemblies are used in their native format, 
and a part’s native geometry engine is accessed through a modeler- 
independent application programming interface (API). In seeking a 
robust and fully automated procedure, the algorithm is based on a 
new physical space manifold triangulation technique specially 
developed to avoid robustness issues associated with poorly condi- 
tioned mappings. In addition, this approach avoids the usual ambi- 
guities associated with floating-point predicate evaluation on 
constructed coordinate geometry in a mapped space. The technique 
is incremental, so that each new site improves the triangulation by 
some well defined quality measure. The algorithm terminates after 
achieving a prespecified measure of mesh quality and produces a tri- 
angulation such that no angle is less than a given angle bound, a, or 
greater than n - 2a. This result also sets bounds on the maximum 
vertex degree, triangle aspect-ratio and maximum stretching rate for 
the triangulation. In addition to the output triangulations for a vari- 
ety of CAD parts, the discussion presents related theoretical results 
which assert the existence of such an angle bound, and demonstrate 
that maximum bounds of between 25° and 30° may be achieved in 
practice. 

1. Introduction 

Mesh generation has long been recognized as a bottleneck in 
the CFD process/ ref ’ ^ The last decade has witnessed a myr- 
iad of international and domestic conferences and sympo- 
siums aimed at focusing research on this impediment. 
Unstructured, hybrid, and Cartesian mesh methods are all 
aimed at simplifying the mesh generation task for complex 
configurations. The success of these approaches is well rep- 
resented in the literature and with an appropriate initial sur- 
face triangulation, the volume mesh generation can generally 
be accomplished in a relatively automated fashion in min- 
utes-to-hours on an engineering workstation. (refs 2_7) As 
faster processors continue to shrink the wall-clock time 
required for both mesh generation and flow solution, the 
man-hour intensive task of extracting an initial surface dis- 
cretization from a CAD geometry promises to become an 
ever larger fraction of the bottleneck. Additionally, if the user 
must be involved in the extraction of surface data from CAD, 
then mesh adaptation - which involves enriching the discreti- 
zation on the body surface - will remain an elusive goal. 


Historically, surface discretization has been one of the least 
automated steps in the numerical simulation cycle, and for 
good reason. Due to its dependence on implicitly defined sur- 
faces and curves, CAD data is by its nature imprecise. Vari- 
ous geometry engines typically demonstrate discrepancies in 
their interpretations of the same entities. As a result, “repair” 
of CAD surfaces has become an area of substantial 
research. 89 This problem is exacerbated when CAD models 
are output in many of the standard formats, since such files 
frequently do not include important topological and construc- 
tion information along with the entity geometry. 

In response to these and other requirements for user assis- 
tance, many in the research and industrial CFD communities 
have adopted an interactive paradigm for surface mesh gener- 
ation. The commercial unstructured mesh generators in refer- 
ences 7,10 and 11 all interact with CAD data through files 
which have been translated from their CAD native environ- 
ment to some standardized format (namely IGES* rcf ” l2 \ 
STEP (ref * 13) or STL (ref - 14) ). 1 

This paper examines an alternative paradigm. The approach 
interfaces with the CAD system in a dynamic manner via 
calls to CAD native routines. By accessing the model in its 
native environment, this approach avoids translation to a for- 
mat which can deplete the model of topological information. 
This is important since it avoids the consistency conflicts that 
can occur when two different geometry engines attempt to 
infer topological information from imprecise data. 

To avoid placing CAD specific calls in the software, we argue 
in favor of wrapping these calls in a standardized Application 
Programming Interface (API) such as CAPRl/ 1 ** 15 ^ This 
library presents a standardized interface to the application 
program for various CAD systems. 2 CAPRI supports a vari- 
ety of operations like truth testing, geometry construction, 
and entity queries. This strategy also divides the task of soft- 
ware maintenance between tasks associated with the CAD 
system and those associated with the surface mesher. 

Maintaining the consistency of the models by direct manipu- 
lation of CAD parts and assemblies is the first step that this 
work takes toward building a robust method for surface trian- 
gulation. A basic premise of this approach is that we resolve 
consistency problems on as simple a model as possible, and 
maintain this consistency as the triangulation evolves and 
becomes more complex. 


1 Recent releases of some of this software now supports “direct” 
interfaces which do read parts in their native formats, however, this 
practice is not the norm. 

2 CAPRI currently supports ProEngineer™ Unigraphics™ and 
SDRC I-DEAS™ with CATIA™ support in beta test. 
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1.1. Abstract Geometric Structures 

Following the approach of Yap, (ref 16) our approach toward 
generating a robust geometric algorithm contends that a geo- 
metric structure, D, consists of four elements: 23 

D = (g, X, <*>(z), z) (1) 

Where the graph, G = (V, E) is a directed set of vertices, V, 
and edges, E.X is a function describing the index labels of 
the graph. O is a geometric operator which represents the 
consistency predicates for the connectivity and is a function 
of the actual coordinates z. G represents a tessellation of the 
vertices and is therefore purely combinatorial. A structure is 
said to be consistent if the predicate O(z) holds. 

As an example, a 2-D Delaunay triangulation algorithm usu- 
ally makes use of an inCircle predicate^ 1 17 \ O inCl>rte which 
establishes G by insisting that the circumcircle of the triangle 
A a b c can contain no other vertex in the graph. If such a pred- 
icate holds for every triangle in G, then this instance of the 
geometric structure D is said to be consistent. This interpreta- 
tion offers direct insights into the formulation of a robust 
algorithm for creating triangulations of CAD volumes. 

1.2. Robustness 

Consistent CAD Models 

The rational B-splines used to describe surfaces in most CAD 
systems are implicitly defined for physical space coordinates 
of the geometry. Therefore, the constructors for vertex geom- 
etry generally require a Newton solve carried to some inter- 
nal tolerance. Since the results of this construction will be 
subject to both tolerance and round-off error, the system may 
then “nudge” the constructed point, z b to some nearby 
exactly representable location (on an integer grid, for exam- 
ple). If the geometry engine’s predicate, ^), for deter- 

mining if Z{ is on a surface, 5 is consistent then it will return 
“true” when later queried if z,- lies on the surface. However, if 
O (Zi S) now represents some user-defined predicate which 
may^e ignorant of the systems construction rules, then it is 
very unlikely to return consistent results. 

The CAPRI API ensures the maintenance of a consistent rep- 
resentation of the model by providing access to a subset of 
the CAD geometry engine’s constructors, queries and predi- 
cates. Our implementation adopts a multi -threaded program- 
ming approach which runs the geometry engine on its own 
thread in order to respond to queries from the main triangula- 
tion thread. 

Physical Space Triangulation 

A variety of existing surface meshing techniques adopt a 
mapped-space approach for generating surface triangula- 
tions. In this approach, a surface and its bounding curves arc 
triangulated in a 2-D parameter space, which may or may not 
have some additional scaling imposed. Reference 18 pro- 


3 Refcrcncel 16] actually writes cq. (1) as D = (G, X, O(z), f) where / 
is a mapping from th^ input parameters c to z, /: 

* c = (cp .... c ) e 9t . 


vides a mathematical description of the Non-Uniform Ratio- 
nal B-Spline (NURBS) surfaces typically used in CAD 
systems. Here we note only that an iterative method is 
required to solve for the physical space coordinates of a posi- 
tion specified on the surface in the parameter space. This pro- 
cess involves division of two (generally) high-order 
polynomials, and is therefore subject to both error associated 
with finite-precision arithmetic and error associated with tol- 
erancing for the convergence of the iterative solve. As a 
result, computed coordinates in the mapped space are neces- 
sarily noisy and cannot be considered exact values. Two con- 
sequences of this approach are: 

1. Since the error bounds on the input, z, are unknown, 
evaluation of the triangulation predicates (e.g. 

O inCircle ) are unlikely to robustly produce consistent 
results (see figure 19 in reference 17; also references 19 
and 20). 

2. The polynomial basis for the NURBS may be high- 
order, and therefore small errors in parameter space may 
produce dramatic results in physical space - even within 
the subspace for which the surface is defined. The likeli- 
hood of encountering poorly conditioned mappings is 
the primary reason that CAD repair software generally 
attempts to renormalize the NURBS surface and recast 
it using basis polynomials with as low an order as possi- 
b l e (re f . 8) 

These two observations motivate an examination of physical 
space triangulation techniques. In this approach, we construct 
a manifold triangulation on the surface, and evaluate the tri- 
angulation predicates in 9t 3 . New sites are constructed by the 
CAD geometry engine and, since this output is consistent 
with the system’s internal predicates, it is considered exact by 
the external predicates of the triangulation algorithm. The 
presentation in §3 emphasizes both minimization and tracing 
of the floating-point error in evaluation of the triangulation 
predicates. Computational requirements for exact evaluation 
are presented. 

2. The CAPRI API 

Our basic approach is to take a crude manifold triangulation 
of each closed volume in the CAD assembly and improve it 
until it either satisfies a preset measure of mesh quality or 
produces a preset number of triangles. A variety of mesh 
quality measures may be defined within this framework, and 
this preliminary investigation examines two such criterion: 

(1 ) the mesh must be free from small angles (sliver triangles); 

(2) edges in the triangulation must not deviate from the 
underlying model by more than a prescribed tolerance. 

2.1. CAPRI Volumes 

CAD entities are accessed through the CAPRI programming 
interfacc. (rcf l5) This API provides a layer of indirection such 
that CAD system specific data may be accessed by an appli- 
cation program using CAD system neutral function calls. 
Figure 1 presents an abstract view of the entities that CAPRI 
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Figure 1- CAPRI data structured demonstrated on a simple volume with a cylindrical cutout. 


provides. A cadjnode is the lowest dimensional entity and 
corresponds to a point in 3-space. A cadjedge has a 
cad_node at both ends. Each edge is directed from its origin, 
O, to its destination, D. Cad_edges are not assumed to be 
simplicial and may follow a general curve in space (see e 5 
and in Fig.l). Each edge is connected to two cad^face 
entities. In general, these faces are composed of several loops 
and are not assumed planar, since they follow the underlying 
parameterization of the surface. A cad_face is composed of 
one or many loops , which are collections of oriented edges. 
Loops are oriented such that the surface of the cad_face lies 
inside them when they are traversed in a counterclockwise 
circuit when viewed from a point outside the solid. This con- 
vention permits holes in a surface to be described by a clock- 
wise loop. In Figure 1, cad_face,/y, consists of loops lj and 
l 2 , each of which is composed of cad_edge entities. The edge 
ordering of l 2 indicates that it is clockwise, and therefore 
describes a hole in fj. Edges and faces have an underlying 
parameterization, and while points may be queried for their 
parametric values ( u , v), details of this ruling are not other- 
wise exposed to the application program. 

2.2. Initial Manifold Triangulation 

A central theme in the present approach is the maintenance of 
a closed volume throughout the procedure. For each volume, 
CAPRI returns a simplicial decomposition of each of the m 
cad_face entities, where {S v S 2 , ...» S m } . Each of 
these triangulations are manifold within their respective 
cad_edges. In addition, an indexing function, X c , is returned 


for each cad_volume. Therefore a simplicial, manifold repre- 
sentation of each cad_volume, Sc may be constructed by tak- 
ing the union of the decompositions of all the cad_face 

entities of a volume, subject to the indexing X^. 

m 

S c = u S, (2) 

i - \ 

Figure 2 displays an example of this initial triangulation for a 
simple part. The manufacturing die shown has 14 cad_face 
entities and the initial triangulation, 5 C , has 270 triangles. As 
is typical, this triangulation is quite irregular, and planar 
regions are decomposed into as few triangles as possible. 
Extremely high aspect ratio triangles are common in these 
boundary triangulations. Figure 2.b labels selected CAD enti- 
ties on this triangulation. Notice that although some cad_face 
sites may be present, this initial triangulation is essentially a 
boundary triangulation and the number of triangles is propor- 
tional to the number of cad_edges. 

Despite the poor quality of the triangulation, the structure in 
Figure 2 has several desirable properties. Namely, it is con- 
sistent, manifold, oriented and closed. We wish to improve 
this triangulation by adding sites on both the cadjedge and 
cad_face entities and by enforcing an external predicate gov- 
erning the type of triangulation. 

3. Mesh Improvement 

Our approach for manifold surface triangulation traces its 
roots to work on quality triangulations of Planar Straight 
Line Graphs (PSLGs) (rcfs * 21, 22 ’ 23) and related work on qual- 




Figure 2. Initial closed, manifold, surface triangulation, S c , of a CAD model for a manufacturing die. Underlying CAD 
entities exposed to application program are labeled in the frame on the right. 
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Figure 3. Constraining edge OD and its diametral circle. The 
edge is “encroached” if any site, p, falls within the diametral cir- 
cle of OD 


ity triangulations of manifold surfaces / ref - 24 ^ Work in this 
field began with the efforts of reference 21 which presented 
an algorithm with both shape and size guarantees. The result- 
ing meshes were size-optimal and had no triangle with an 
aspect ratio greater than 5. In this context, the aspect ratio 
AR, of a triangle, is defined as the length of the longest edge 
divided by that of the shortest one. One can show that if a is 
the smallest angle of a triangle, then 


1 


< AR <i 


(3) 


| sinoc| | sina| 

Therefore a is frequently used to describe the quality of a 
given triangle. 


Before presenting the manifold triangulation technique, this 
section first recounts a related algorithm for quality triangula- 
tion for PSLGs from reference 23. It then presents a funda- 
mentally similar algorithm for triangulating curved surfaces 
and notes which aspects of the PLSG method have been 
relaxed in the extension. 


3.1. Quality Triangulation of PSLGs 

While the manifold surface triangulation technique of Rup- 
pert (rcf - and the PSLG method of Chew (rcf ‘ 24) are similar 
in many respects, our manifold technique follows Ruppert’s 
approach more closely. Section 3.2 addresses some of the 
reasoning behind this choice. 


predicate to an edge, e. The (~) superscript reminds us that 
since this predicate is part of our triangulation algorithm, it is 
not native to the CAD system. 

The presentation of reference 23 recovers the constraining 
edges of the triangulation as the algorithm advances. To clar- 
ify its relation with our manifold triangulation technique, we 
recast the original algorithm assuming that it begins with a 
constrained boundary triangulation of the input vertices, V, of 
the constraint edges. Furthermore, this initial triangulation is 
assumed to be the constrained Delaunay triangulation of the 
input sites, C2)%V). 

The algorithm is quite elegant in that it consists of only two 
major operations: 

1. Split a constrained edge: Add a site near 5 the mid-point 
of a constraining edge, and replace the original edge 
with the two new edges in the constraint list. 

2. Split a triangle: Add a site to the circumcenter of a trian- 
gle, t. Note that if the triangle is obtuse, this site will not 
fall within t. 

Algorithm Q: Quality triangulation of a PSLG 

Input: Planar Straight Line Graph, X , with input vertices, V in . 

Target angle, a. 

Output: ClXI{y out ) with all angles > a. 

Initialize: Compute CTH{V in ). Build minimum angle priority 
queue, PQ m i n with t P Q denoting triangle at head of 
queue, having min angle Q P q 

1. Apply Encroached^) to all constraint edges: 

While(any constraining edge is encroached)! 

Split constrained edge. Update CZXZ(V), Update PQ m i n . 

1 

2. While (0po < °0{ 

2.a Let p be the circumcenter of t P Q. 

2.b If ip encroaches any constraining edge, e) 

2.c Split constrained edge. Update Cin{ V), Update PQ m i n . 

2. d Else Split triangle t P Q. 

Add p to V. Update UZ?7(V), Update PQ m i n . 

} 

3. Output CD%V). 


An essential feature of the algorithm is the notion of an 
encroached constraining edge. As illustrated in Figure 3, a 
constraining edge, OD , is said to be encroached upon if any 
other site (visible to OD ) lies within the diametral circle of 
the edge. 

If one recalls that the circumcenter of a right triangle falls on 
the hypotenuse, then its easy to show that for a triangulation 
which is Delaunay or locally maxmin, a predicate for 
encroachment may be formulated as a vector dot product. 
Thus for similarly sized 4 floating-point data with p-bit 
significands, this predicate can be evaluated exactly in a reg- 
ister with a 2p-bit significand (rcf 20) . In a practical sense, this 
implies that as long as the edges are small by comparison to 
their distance from the origin, this predicate will be exact if 
computed in double-precision, using single-precision data. 

In the algorithms, ^ Encroached^) denotes application of this 


While simple, reference 23 proves that Alg. Q produces tri- 
angulations with the following desirable properties. 

1. Quality: An angle bound 20.7° is guaranteed, and values as 
high as 30° may be achieved in practice. 

2. Output Size: The size of the output triangulation is within 
a constant of the optimal number of triangles required to sat- 
isfy the angle criteria. 

3. Size Optimality: Small input constraint edges are sur- 
rounded by proportionally small triangles. Nearby triangles 
have similar sizes, and the size variation of triangles in the 
mesh is proportional to the distance between them. 

3.2. Difficulties on Curved Surfaces 

In reference 24, Chew presents a quality triangulation tech- 
nique which is closely related to Alg. Q. This work raises a 
number of difficulties associated with the extension of the 


4 The qualifier “similarly sized” is necessary to guard against the 

case where a coordinates of one point is less than half that of 

another. Extended precision would be required in such a case. 5 We use the modified insertion strategy presented in reference 23. 
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Figure 4. An example of nonintuitive consequences of a straight- 
forward interpretation of the inCircle PSLG predicate on curved 
surfaces. With distances measured using the geodesic distance 
along the surface, the loci of points equidistant from p reaches 
around the spire but does not include its tip. 

PSLG method to curved surfaces, alg. 


tive method and will therefore be computationally intensive. 

(3) Once this intersection is successfully located on the sur- 
face, one must determine which triangle the point falls inside. 
Since the triangulation only matches the surface at the verti- 
ces, ambiguous situations may arise when two triangles 
claim ownership of the same site (near a ridge, for example). 

(4) On a curved surface, the circumcenters of two triangles 
with a shared edge may not be consistent. Specifically, when 
6 inCircle tests one °f the triangles against the opposite ver- 
tex of the other, the results may not agree when the roles of 
the triangles are reversed. Lemma 5 in reference 24 shows 
this situation will arise if the normal vectors of the triangles 
vary by more than rc/2 . A successful algorithm must guard 
against O inCircle becoming inconsistent. 

3.3. A Physical-Space Surface Meshing Algorithm 

These outstanding ambiguities and computational expense 
motivated a search for a more manageable algorithm. Our 
method makes two fundamental changes. 


Alg. Q has two salient aspects. (1) The triangulation is con- 
strained Delaunay. (2) New sites are added at circumcenters. 
Chew observes that a straightforward definition of a circle on 
a surface is the loci of points on the surface which are equi- 
distant from another point on the surface, where all distances 
are measured using the geodesic distance along the surface. 
While straightforward, this definition is problematic. Dis- 
tances along the surface must be measured in physical space, 
and will therefore be expensive to compute on NURBS sur- 
faces. In addition, due to the inherent error in finding the 
coordinates of a point on such a surface, robust predicates 
based upon this definition will be difficult to formulate. 
Finally, Chew notes that this definition has less subtle and 
non-intuitive consequences. As shown in Figure 4, a circle 
whose center lies near the base of a sharp spire, for example, 
may reach completely around the spire without also includ- 
ing the tip of the spire. 

To circumvent such difficulties, the method in reference 24 
makes use of an alternative definition of a circle. In the plane, 
the three vertices of a triangle define a unique circle. In three 
dimensions, however, an infinite family of spheres may be 
passed through those three points. Connecting the line 
through the centers of this family of spheres and intersecting 
this line with the surface identifies a particular sphere in this 
family. The circumcenter of the triangle may be defined to be 
the loci of points at the intersection of this particular sphere 
and the surface. Once this sphere is found, then the 
O inCircle triangulation predicate may be evaluated by sim- 
ply computing distances in three dimensions as a vector mag- 
nitude. 

Despite the effort, problems still exist with this predicate. (1) 
If the triangulation does not resolve the underlying surface 
closely enough, the line of circumsphere centers will not nec- 
essarily intersect the surface. Alternatively, it may also inter- 
sect the surface in multiple locations. (2) Computing the 
intersection of this line with the surface will require an itcra- 


• To avoid the ambiguity associated with the definition of 
6 inCircle i we d° not attempt a Delaunay triangulation of 
the curved surface. Instead we seek a triangulation which is 
everywhere locally maxmin. While this may seem to consti- 
tute a dramatic relaxation, recall that the goal is a practical 
algorithm, and we have no reason to prefer strictly Delaunay 
output triangulations. In addition, since one property of a 
Delaunay triangulation is that it is maxmin, this choice is 
worth investigating. 

• When all angles of a triangle are acute, the circumcenter 
falls within the triangle. Ownership of the new site can then 
be uniquely assigned to this triangle. However, when a trian- 
gle is obtuse, ownership can become less clear. Therefore we 
make a simple choice: When a triangle is obtuse, we insert 
the new site at the centroid of the other triangle sharing the 
edge opposite the obtuse angle. Figure 5. illustrates this 
insertion rule. If tj is an obtuse triangle, then the new sites are 
added at the centroid of t opp . 

While admittedly ad hoc, the modified insertion strategy is 
not as arbitrary as it may initially seem. As a particular angle 
of ty opens up in the transition from acute to obtuse, the cir- 
cumcenter will pass from within t t to within t opp . The cen- 
troid of t opp may therefore be thought of as an approximate 




tj is acute /, is obtuse 

Figure 5. Sites arc added at the centroid of a triangle if the trian- 
gle is acute. Otherwise they are added at the centroid of the tri- 
angle opposite the obtuse angle. 
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circumcenter of Denoting the radius of the circumcircle of 
t t as R, it can be shown that if one chooses an approximate 
circumcenter within a distance fR of the true circumcenter, 
then an angle bound of 

a< asin(— 2 ^-) (4) 

may still be achieved. Reference 24 contends that even using 
approximate circumcenters, angle bounds of 30° may be 
achieved in practice. 

With these changes, the manifold triangulation algorithm 
becomes: 


incremental insertion strategy of reference 25. However, 
since the triangulation is no longer Delaunay, both forward 
and reverse propagation is necessary after edge swaps. Sites 
from edge or triangle splitting must be projected to their 
actual locations on the underlying surface, and the native 
constructors are used for this through the CAD API. 
Reference25 gives a modified point placement strategy for 
splitting encroached edges. Our implementation of Alg. M 
adopts this strategy without modification. 

3.4. Triangulation Predicates 


Algorithm M: 

Quality Manifold triangulation of a CAD volume. 

Input: Underlying CAD volume, P, with initial triangulation 
S c , and input vertices, V in . 

Target angle, a. 

Output: CX$&V out ) with all angles > a. 

Initialize: Compute constrained locally maxmin triangulation 
CX%l[V in ), using cad_edge entities as constraints. 
Build minimum angle priority queue, PQ m (n with t P Q 
denoting triangle at head of queue, having min angle 
~ ®PQ- 

1. Apply O Encroached to all constraint edges: 

While(any constraining edge is encroached)! 

Split constrained edge. Update CX9&V), Update PQ min . 

2. While (Q PQ < a){ 

2.a If ( t P Q is not obtuse) 

Assign t := t P Q, with circumcenter p. 

Else Assign t:= t opp with centroid p. 

2.b If ip encroaches any constraining edge, e) 

2.c Split constrained edge. Update CX9&V), Update PQ min . 

2. d Else Split triangle t\ 

Add p to K Update CtiAtfV), Update PQ min . 

1 

3. Output CT^V). 

When the algorithm terminates, all angles are greater than a. 
Thus, it recovers the properties of quality and size-optimality 
cited after the presentation of Alg. Q in §3.1. The bound on 
output size, however, depends on the site insertion strategy 
always inserting new vertices within the circumcircle of 
t P Q , (ref - 23) and the modified strategy will not always guaran- 
tee this, thus the algorithm sacrifices strict proof of this prop- 
erty. 

The new algorithm requires initialization with a transforma- 
tion of the part’s original constrained manifold triangulation 
S c to one which is everywhere locally maxmin. If we assume 
the existence of a triangulation predicate which can be 
enforced for every pair of triangles sharing an edge in the tri- 
angulation (similar to the application of Q>inCircle )> then 
this initialization may be performed with edge sweeps fol- 
lowed by edge-swapping when a violation is encountered. 
While the possibility of multiple sweeps makes this a seem- 
ingly inefficient approach, we recall that the initial complex- 
ity of S c is only proportional to the number of cad_cdgcs in 
the geometric structure. Thus this simplistic approach is not a 
problem. 

Site insertion proceeds in a manner comparable to the PSLG 
algorithm. A convenient implementation makes use of the 


Algorithm M rests on two new triangulation predicates. The 
first tests a triangle for an obtuse angle, 0^(0 and the sec- 
ond, i>xN(e) , tests if the edge, e, shared by any two triangles 
maximizes the minimum angle in both triangles. This is 
accomplished by comparison with a swapped edge, e ', which 
connects the opposite vertices of the two triangles. Our 
approach hinges on the hope of evaluating these predicates 
robustly in physical space. 

Both of these predicates can be formulated with direct mea- 
surement of the angles in a mesh. Since a triangle has three 
points, it uniquely defines a plane in three dimensions. Angle 
measurement in a plane is unambiguous - despite the fact 
that the triangles form a discrete manifold which is of lower 
dimension than the surrounding space. Numerically, we 
recall that since all vertex locations returned by the CAD 
engine are considered exact, the error bound on the input is 
identically zero. 

i>ob ( 0 ft applied by a logical or accumulation of 
4> zj-\,jj+ 1 , where j is a cyclic index running over the verti- 
ces of f, j e { 1, 2, 3} . This angle predicate is formed by com- 
parison of the square of the edge opposite j to that of a 
hypothetical edge formed by assuming that the edges of t 
incident on j form a right angle in the plane of t. 

^ob^m) = Kb( Zl2i ) V<t> 0 fcU23l) V<j> oi (Z312) (5) 

where 


T if b-^-i| 2 + K.>- v J 2 >K*>-V; -if (6) 

F otherwise 

Notice that eq. (6) requires only subtraction and multiplica- 
tion of data which is known exactly. Thus the computational 


*4 



Figure 6. Construction for \d\ 2 , the square of the magnitude of 
the difference of unit vectors incident on vertex j. 
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requirements for exact evaluation are the same as for evalua- 
tion of the dot product for <&Encroached(e) in §3.1. 

Evaluation of 6 XAr (<?) requires a more direct method of angle 
measurement. Figure 6 shows a construction for this mea- 
surement. Recalling that sin( ) is monotone over the interval 
[0, rc/2] the construction in the figure shows that the differ- 
ence of the unit vectors of the edges incident upon any ver- 
tex, j, is sufficient to define a vector 3 whose magnitude 
varies monotonically with the angle formed by the edges 
incident upon j. As in the preceding two predicates, computa- 
tionally, it is sufficient to evaluate only the square of this 
magnitude. 

The computational requirements for this predicate are some- 
what less simple. The computation of unit vectors requires 
robust computation of the inverse of a vector magnitude 
l/|v| . Thus, exact evaluation of <&xN(e) requires software 
arithmetic, like the packages in references 26, 27 or 28. 
Although this need must be viewed as a drawback, we note 
that it is confined to a single predicate, and to a relatively 
simple expression. Moreover, since the input geometry is 
exact, exact computation remains a feasible strategy. In the 
preliminary results shown in §4, all computation was per- 
formed using only double-precision floating- point hardware, 
and the option of software arithmetic was not pursued. 

3.5. Edge refinement 

Algorithm M drives small angles out of the evolving triangu- 
lation. In addition we wish to satisfy some edge-based crite- 
rion like a chord-height tolerance. After initially creating a 
triangulation free of small angles we apply an edge-refine- 
ment procedure to enforce such requirements. Algorithm E 
considers a generic edge-based scalar y(e). In our implemen- 
tation chord-height is defined as the square of the distance 
from the middle of an edge to the corresponding location on 
the actual surface of the model (provided through the CAD 
API by the geometry engine). 

Algorithm E: Edge refinement of Manifold Triangulation. 

Input: Underlying CAD volume, P, with current triangulation, 
CX9\&V), and vertex set, V. 

Edge criteria y. 

Output: CX?4y out ) with all edges satisfying y(e) < y. 

Initialize: Build priority queue, PQy with epQ denoting edge at 
head of queue, having y (e P Q). 

1 . Apply O Encroached to all constraint edges: 

While(any constraining edge is encroached)! 

Split constrained edge. Update CXi A&V), Update PQ y . 

} T 

2. While (y(epQ) > y){ 

2. a Let p be midpoint of e. 

2.b If {p encroaches any constraining edge, e) 

2.c Split constrained edge. Update OC9&V), Update PQy. 

2. d Else Split edge e : 

Add p to e. Update CX9&V), Update PQ r 

} 

3. Output CXOiY). 

Assuming that y(e) is a static criterion, like a chord-height 
tolerance, Alg. M can then be re-applied to CX9&V) to 
remove any small angles created during edge refinement and 


swapping. 

4. Results and Discussion 

This section presents example meshes on several CAD parts 
from a variety of sources. All the example parts were read in 
their native CAD file format using the CAPRI API without 
special treatment. The investigations focus on examination of 
issues raised in the presentation of the triangulation algo- 
rithm in §3.3 and the edge refinement strategy from §3.5. 

4.1. Minimum Angle Bound 

In §3.3, Algorithm M was presented without firm proof of 
termination. Moreover, the discussion noted that the modified 
site insertion strategy for obtuse triangles violates one of the 
assumptions that establishes a bound on the output size of the 
mesh in the PSLG method. It is therefore necessary to dem- 
onstrate the performance of Alg. M to show that it both ter- 
minates and produces meshes with an economy similar to 
that of the PSLG method upon which it is based. 

Figure 7 contains a histogram of the evolution of the smallest 
angle in the mesh as Alg. M proceeds on the manufacturing 
die example problem used in earlier illustrations. While this 
curve is far from monotone, it clearly displays the steady 
improvement of the minimum angle in the mesh as the algo- 
rithm proceeds. The steep initial rise indicates rapid annihila- 
tion of extremely small angles in the mesh, and the mesh 
achieves a minimum angle of almost 29° by the end of the 
histogram. 

The dashed line at 25° highlights the first time that all angles 
in the mesh exceeded this value. Tracing this value on the 
abscissa shows that setting the angle bound, a, to 25° will 
cause Alg. M to terminate after generating 2606 triangles. 
Figure 8 shows the resulting triangulation. As discussed in 



0 — * — J 

0 2000 4000 6000 8000 10000 

Number of Triangles in Manifold 

Figure 7. Histogram of minimum angle during mesh evolution 
using Alg. M (without the edge refinement of §3.5) for the man- 
ufacturing die example presented earlier. In this example, a mesh 
with a minimum angle of 25° would contain 2606 triangles. 
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Figure 8. Quality manifold triangulation for manufacturing die 
example generated with Alg. M. in §3.3. Mesh improvement 
terminated after generating 2606 triangles when the mini- 
mum angle in the triangulation reached 25°. Chord-height 

§3.1 and §3.3, the presence of an angle bound ensures that 
small features are surrounded by proportionally small trian- 
gles (see the inset frames in figure 8), and that the mesh 
length-scale varies smoothly over the part. The smallest 
angle in the mesh is 25.02°, which corresponds to a maxi- 
mum aspect ratio of between 2.36 and 4.7 (by eq. (3)). 

While the histogram in Figure 7 shows a steady increase in 
the minimum angle, there is an irregular array of downward 
spikes in the profile. In the presentation of the PSLG algo- 
rithm, reference 23 included a similar histogram and noted 
the same characteristic. Consider the two triangles t x and t 2 
shown at the left in Figure 9, Assume that Q x is the smallest 
angle in the mesh lying in triangle t } . Furthermore, assume 
that s face neighbor, t 2 , has a small angle 0 2 opposite the 
shared edge which is only slightly larger than Q x . Site p will 
get added at 7/s circumcenter, improving d x to 2 0^ After 
application of the maxmin predicate <b X N on the shared edge, 
the swapped configuration at the right of Fig.9 may occur. 
This configuration includes two new triangles with minimum 
angles 0 3 and 0 4 either of which may now be the smallest 
angle in the mesh and may actually be smaller than the origi- 
nal angle 0j. 

With this behavior understood, Figures 10 and 11 present 
angle histograms and example meshes for a more compli- 



Figure 9. Mechanism responsible for downward spikes in histo- 
gram of minimum angle as Alg. M proceeds. If 0] is initially 
the smallest angle in the mesh, angles 0^ and 0 4 may be smaller 
after insertion of p and enforcement of <P xyv by edge-swapping. 


cated assembly of parts. Alg. M was run on CAD parts for 
the main element of a transport wing, and a flap element for 
the same wing. The main element consisted of 224 rational 
B-spline curves and 36 trimmed NURBS surfaces. The flap 
contained 31 rational B-spline curves and 10 trimmed 
NURBS surfaces. 

The crosshairs on the curves in figure 10 show that with a 25° 
angle bound, Alg. M will produce 20846 triangles on the 
main element and 15334 triangles on the flap. Figure 11 dis- 
plays these triangulations. 

The histograms in Figure 10 bear close resemblance to the 
one presented for the simple die example shown in figure 7. 
All of these profiles are characterized by a sharp initial angle 
improvement and then a rolling-off as the minimum angle 
climbs above about 20°. All three profiles exceed 27°, but 
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Figure 10. Angle histograms for Alg. M on the main element of a 
transport wing(uppcr), and on the main flap (lower). An angle 
bound of 25° would produce triangulations with 20846 and 
15334 triangles on the main element and flap respectively. 
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Figure 11. Bounded angle triangulations of main element of a 
transport aircraft wing and flap generated by Alg. M in §3.3. 
Minimum angle 25°, 20846 triangles on wing, 15334 triangles 
on flap. 

reaching 30° seems unlikely. 

The abscissa values in figure 10 indicate relatively large tri- 
angulations as compared with the die example. The size-opti- 
mality property cited in §3.1 and §3.3 implies that triangle 
size must vary smoothly between different sized constraints. 
In both of these examples the smallest constraint in the 
CAPRI triangulation was a factor of 10 4 smaller than the 
largest constraining edge. In addition, practical experience 
with the algorithm indicates that larger initial triangulations, 
S c , (eq. (2)) have a proportionally slower initial rise in their 
angle profiles. This seems reasonable since if the CAPRI tri- 
angulation is complex there may initially be many bad angles 
which need improvement. Often CAD parts are unnecessarily 
complex for reasons dating to the specific events in their cre- 
ation. Our experience with CAD repair software indicates 
that CAPRI generally produces less complex (sometimes by 
an order of magnitude) initial triangulations if the parts have 
been processed by CAD repair tools. Since tolerance to poor 
CAD parts is one type of robustness that we seek in this 
research, neither part shown in Fig.l 1 underwent such repair 
prior to creation of these triangulations. 

With the minimum mesh angle restricted to 25°, triangles are 
again restricted to aspect ratios less than 4.7. Such an isotro- 
pic mesh is very inefficient at meshing features with curva- 
ture in only one dimension. Therefore, if one intends to 
produce meshes for viscous computation, a substantially 
smaller angle bound may be appropriate. 

4.2. Chord-Height Refinement 

Enforcement of the angle bound does not guarantee that the 
edges are refined when they are sufficiently far from the 
underlying surface. Alg. E in §3.5 presented an edge refine- 
ment strategy for automatically breaking edges whose mid- 
points are far from the geometry. Figure 12 shows an interior 
comer of the die example after the imposition of a chord- 
height tolerance of 10~ 4 times L, where L is the maximum 



Figure 12. Interior comer of manufacturing die example trian- 
gulation after imposition of a chord-height tolerance of 
l.xlO -4 normalized by the maximum outer dimension of the 
part. Edge refinement was performed using Alg. E. in §3.5 

outer dimension of the part. Away from the curved interior 
comers, the triangulation remains unchanged since the faces 
of the die are planar. 

5. Conclusions and Future Work 

5.1. Conclusions 

This paper examined the direct use of CAD geometry in sur- 
face meshing. It argued in favor of accessing CAD parts and 
assemblies in their native format, without translation and for 
the use of CAD-native predicates and constructors in mesh 
generation. The discussion also contended that since physi- 
cal-space mesh generation techniques permit exact numerical 
computation, they can be more robust than their mapped- 
space counterparts. 

The feasibility of this approach was demonstrated by applica- 
tion to a novel, physical-space, curved surface meshing algo- 
rithm to several CAD parts. Model data was accessed in its 
native format using the CAPRI API. This triangulation algo- 
rithm maintained a locally maxmin triangulation and used an 
approximate circumcenter site insertion strategy. All triangu- 
lation predicates were formulated for evaluation in physical- 
space, and examples demonstrated that the method produced 
a bounded aspect ratio manifold triangulation. The computa- 
tional requirements for exact evaluation of all triangulation 
predicates were discussed. The algorithm was demonstrated 
on CAD parts of varying complexity and reliably produced 
triangulations with minimum angles in excess of 27°. 

5.2. Future Research 

• Although Alg. M. produces meshes with generally smooth 
length-scale variation, there is occasionally a discernible 
irregularity in the triangulations. This behavior becomes 
more pronounced when higher angle bounds are specified. 
Similar behavior has been noted in the PSLG algorithm 23, 
although in two dimensions, the behavior seems less pro- 
nounced. One possible source for this is the abrupt change in 
insertion location from circumcenter to centroid (of t opp ) if 
an obtuse triangle is encountered. Alternative strategies 
should be investigated since this behavior can degrade the 
efficiency of the triangulations. 

• The initial triangulation returned by CAPRI depends 
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strongly upon the “history” of the CAD part. Such effects are 
evident in the triangulation of the main wing element in fig- 
ure 11, for example, where the outboard portion of the flap 
cut-out is excessively resolved due to relics of the part’s cre- 
ation process. The only effective strategy we’ve found for 
avoiding this is to “clean” the geometry using CAD repair 
software. Alternative approaches should be investigated. 

6. References 

1 Cosner, R.: Issues in Aerospace Applications of CFD Analysis. 

AIAA Paper 94-0464, Jan. 1994. 

2 Aftosmis, M.J.; Berger, M.J.; Melton, J.E.: Robust and Effi- 

cient Cartesian Mesh Generation for Component-Based 
Geometry.” AIAA Paper 97-0196 , Jan. 1997. 

3 Lohner.R.: Finite Element Methods in CFD: Grid Generation, 

Adaptivity and Parallelization. AGARD Special course on 
unstructured Grid methods for Advection dominated Flows, 
AGARD-R-787, May, 1992. 

4 Marcum, D. L.; and Weatherhill, N. P.: Unstructured Grid Gen- 

eration Using Iterative Point Insertion and Local Reconstruc- 
tion.” AIAA J. vol. 33, no, 9, Sept. 1995. 

5 Pirzadeh, S.: Unstructured Viscous Grid Generation by the 

Advancing Layers Method. AIAA Paper 93-3453-CP. Jun. 
1993. 

6 Wang, Z. J.; Przekas, A.; and Hufford, G.: Adaptive Cartesian/ 

Adaptive Prism Grid Generation for Complex Geometry. 
AIAA Paper 97-0860 , Jan. 1997. 

7 ICEM CFD/CAE V.3.1.3: User Manual, Volume 2, Unstruc- 

tured Grid Generation. ICEM Systems Inc. Berkeley CA, 
Nov. 1994. 

8 Bohn J. W.; and Zozny, M. J.: Automatic CAD-model repair: 

Shell -closure. Proc. Symp. On Freeform Fabrication. Dept, 
of Mech. Eng.; Univ. of Texas at Austin, 1992. 

9 Gu£ziec, A.; Taubin, G.; Lazarus, F.; and Horn, W.: Cutting 

and Stitching: Efficient Conversion of a Non-manifold 
Polygonal Surface to a Manifold”, IBM-RC-20935. IBM 
Research Division, Yorktown Heights, NY, Jul. 1997. 

10 VGRID: An Unstructured Tetrahedral Grid Generator Based 

on the Advancing Front Method. Vigyan Inc. Hampton VA. 
http://www.vigyan.com, 1998. 

1 1 CFD-GEOM User’s Manual. CFD Research Corp., Huntsville 

Al, 1997. 

12 Reed, K.: The Initial Graphics Exchange Specification (IGES) 

Version 5.1. Sept. 1991. 

1 3 Industrial Automation Systems and Integration - Product Data 

Representation and Exchange - Part 1 :Overview and Funda- 
mental Principles. ISO/TR 10303-L International Standards 
Org. Geneve, Switzerland, 1994. 

14 3D Systems Inc.: Stereolithography Interface Format Specifi- 

cation. 1988. 

15 Haimes, R.; and Follen, G.: Computational Analysis PRogram- 

ming Interface. Proceedings of the 6th International Confer- 
ence on Numerical Grid Generation in Computational Field 
Simulations , Eds. Cross, Eiseman, Hauser, Soni and Thomp- 
son, July 1998. 

16 Yap, C.: Robust Geometric Computation. CRC Handbook of 

Discrete and Computational Geometry Eds. Goodman J.E.; 
and O’Rourke, J., CRC Press, Boca Raton FI.; 1997. 

17 Shewchuk, J.R.: Adaptive Precision Floating-Point Arithmetic 

and Fast Robust Geometric Predicates. CMU-CS-96-140. 
Carnegie Mellon Univ. School of Computer Science, May, 
1996. 

18 Farm, G.: Curves and Surfaces for Computer Aided Geometry 

Design, A Practical Guide. Academic Press, 1990. 

19 Michclucci, D.: The Robustness Issue. Internal Report, Labo- 

ratoirc d’ Image dc Synlh6se de St Etienne, France. Sec hup:/ 
/www.cmsc.fr/~michcluc/cnglish/michclucci.html 


20 Goldberg, D.: What Every Computer Scientist Should Know 

About Floating-Point Arithmetic. ACM Comput. Surveys, 
vol 23, no. 1, pp. 5-48, Mar. 1991. 

21 Bern, M.; Eppstein, D.; and Gilbert, J.R.: Provably Good Mesh 

Generation, in Proceedings of the 31st Annual Symposium on 
Foundations of Computer Science . pp. 231-241, IEEE, St. 
Louis, Missouri, Oct. 1990. 

22 Chew, L.P.: Guaranteed-Quality Triangular Meshes. Technical 

Report, TR-89-983. Computer Science Dept, Cornell Univ. 
Apr. 1989. 

23 Ruppert, J.: A Delaunay Refinement Algorithm for Quality 2- 

Dimensional Mesh Generation. J. of Algorithms, vol. 18, no. 
3, pp. 548-585, May 1995. 

24 Chew, L.P.; “Guaranteed-quality mesh generation for curved 

surfaces.” Proc. of the Ninth Annual Symposium on Compu- 
tational Geometry , 274-280. ACM, 1993. 

25 Green, P.J.; and Sibson, R.: Computing the Dirichlet Tessella- 

tion in the Plane. The Computer J. vol. 2, no. 21, pp. 168- 
173, 1977. 

26 Real/Expr homepage: Source Code, Documentation, Examples 

and Literature, http://simulation.nyu.edu/projects/exact, 
1996. 

27 Ouchi, K.: Real/Expr: Implementation of an exact computation 

package. Master’s Thesis, New York Univ., Dept, of Comp. 
Sci., Courant Inst. NY. Jan. 1997. 

28 Trimaran homepage: http://react-ilp.nyu.edu. 1998. 


10 



REPORT DOCUMENTATION PAGE 


Form Approved 
QMS No 0704-0188 


Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time tor reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 
collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for information Operations and Reports, 1215 Jefferson 
Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-Q188), Washington, DC 20503. 


1. AGENCY USE ONLY (Leave blank) 


2. REPORT DATE 

August 1999 


REPORT TYPE AND DATES COVERED 

Technical Memorandum 


4. TITLE AND SUBTITLE 

On the Use of CAD-Native Predicates and Geometry in Surface Meshing 


5. FUNDING NUMBERS 


6. AUTHOR(S) 

M. J. Aftosmis 


522-31-12 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

Ames Research Center 
Moffett Field, CA 94035-1000 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 


A-99V0022 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


NAS A/TM-1 999-208782 


11. SUPPLEMENTARY NOTES 

Point of Contact: M. J. Aftosmis, Ames Research Center, MS T27B-2, Moffett Field, CA 94035-1000 
(650) 604-4499 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified — Unlimited 
Subject Category 61 

Availability: NASACASI (301) 621-0390 


12b. DISTRIBUTION CODE 


Distribution: Standard 


13. ABSTRACT (Maximum 200 words) 

Several paradigms for accessing computer-aided design (CAD) geometry during surface meshing for 
computational fluid dynamics are discussed. File translation, inconsistent geometry engines, and normative 
point construction are all identified as sources of nonrobustness. The paper argues in favor of accessing 
CAD parts and assemblies in their native format, without translation, and for the use of CAD-native predi- 
cates and constructors in surface mesh generation. The discussion also emphasizes the importance of 
examining the computational requirements for exact evaluation of triangulation predicates during surface 
meshing. 


14. SUBJECT TERMS 



Triangulation algorithm, CAD, Curved surface 


17. SECURITY CLASSIFICATION 

18. SECURITY CLASSIFICATION 

19. SECURITY CLASSIFICATION 

OF REPORT 

OF THIS PAGE 

OF ABSTRACT 

Unclassified 

Unclassified 



NSN 7540-01 -280-5500 


15. NUMBER OF PAGES 

15 

16. PRICE CODE 

A03 


Standard Form 298 (Rev. 2-89) 

Prescribed by ANSI Std. Z39-10 
290-102 









